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Abstract 

The wetlands of the Qinghai-Tibetan Plateau are believed to play an important role in global nutrient cycling, but the 
composition and diversity of microorganisms in this ecosystem are poorly characterized. An understanding of the effects of 
geography and microtopography on microbial populations will provide clues to the underlying mechanisms that structure 
microbial communities. In this study, we used pyrosequencing-based analysis of 16S rRNA gene sequences to assess and 
compare the composition of soil microbial communities present in hummock and hollow soils from three wetlands 
(Dangxiong, Hongyuan and Maduo) on the Qinghai-Tibetan Plateau, the world's highest plateau. A total of 36 bacterial 
phyla were detected. Proteobacteria (34.5% average relative abundance), Actinobacteria (17.3%) and Bacteroidetes (11%) 
had the highest relative abundances across all sites. Chloroflexi, Acidobacteria, Verrucomicrobia, Firmicutes, and 
Planctomycetes were also relatively abundant (1-10%). In addition, archaeal sequences belonging to Euryarchaea, 
Crenarchaea and Thaumarchaea were detected. Alphaproteobacteria sequences, especially of the order Rhodospirillales, 
were significantly more abundant in Maduo than Hongyuan and Dangxiong wetlands. Compared with Hongyuan soils, 
Dangxiong and Maduo had significantly higher relative abundances of Gammaproteobacteria sequences (mainly order 
Xanthomonadales). Hongyuan wetland had a relatively high abundance of methanogens (mainly genera Methanobacterium, 
Methanosarcina and Methanosaeta) and methanotrophs (mainly Methylocystis) compared with the other two wetlands. 
Principal coordinate analysis (PCoA) indicated that the microbial community structure differed between locations and 
microtopographies and canonical correspondence analysis indicated an association between microbial community 
structure and soil properties or geography. These insights into the microbial community structure and the main controlling 
factors in wetlands of the Qinghai-Tibetan Plateau provide a valuable background for further studies on biogeochemical 
processes in this distinct ecosystem. 
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Introduction 

Soil microbial communities in wetland systems play an 
important role in biogeochemical cycles and are crucial to the 
functions of wetland systems [1]. In recent years, owing to the 
advantages of next-generation sequencing, research on microbial 
diversity in various natural wetlands has rapidly developed. For 
example, the 16S rRNA tag-encoded pyrosequencing approach 
has been adopted to study microbial diversity in different natural 
wetlands, such as a Sphagnum-dominated northern wetland of 
Russia [2], marshes in eastern Massachusetts [3], peatlands in the 
Glacial Lake Agassiz region [4], and freshwater wetlands in 
northern Virginia [5]. Although considerable progress has been 
achieved in characterizing the microbial diversity in natural 
wetlands, further research is still needed. 

The Qinghai-Tibetan Plateau (av. 4000 m a.s.l.) is the largest 
and highest plateau on Earth. The area of natural wetlands on the 



Qinghai-Tibetan Plateau is estimated at 13.3xl0 4 km 2 . In contrast 
to boreal wetlands located at high latitude regions, wetlands in this 
region are at low latitude and high altitude. Additionally, they 
have distinguishing features of hummock-hollow microtopography 
and Sedge-dominated vegetation [6]. Wedands on the Qinghai- 
Tibetan Plateau are very sensitive to climate change and they play 
a vital role in terrestrial carbon storage [7]. Also, these wedands 
provide many important ecological services, such as ensuring the 
productivity of rangelands and providing grazing ground for 
thousands of livestock animals. Despite the importance of these 
wetlands, the composition of soil microbial communities in this 
distinct ecosystem has not been fully investigated. The microbi- 
ological research in this region has been primarily restricted to 
specific functional groups of microorganisms, such as methanogens 
[8,9] and methanotrophs [10,11], or microbes involved in other 
specialized processes [12,13]. In addition, a few cultivation- 
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independent studies using denaturing gradient gel electrophoresis 
(DGGE) fingerprinting and limited clone library construction have 
been carried-out to assess the bacterial diversity in wetlands of this 
region [14]. 

Due to practical limitations of the techniques previously 
employed, limited phylogenetic resolution and coverage of the 
microbial diversity in wetland environments on the Qinghai- 
Tibetan Plateau is available. Furthermore, many important 
aspects are still not investigated. First, despite numerous wetlands 
unevenly distributed on the Qinghai-Tibetan Plateau [15], the 
differences in microbial communities between these wetlands at 
different geographical locations are mostly unexplored. Second, no 
study has focused on variations in the microbial community 
structures between hummock and hollow soils, which is the 
common microtopography in wetlands on this plateau [6] . Lastly, 
recent studies using molecular techniques have provided evidence 
that environmental factors, particularly soil properties, influence 
the bacterial community structure in several terrestrial and aquatic 
systems including wetland environments [16-19]; however, this 
has not been investigated in wedands of the Qinghai-Tibetan 
Plateau. 

The aim of the present study was to detect the composition and 
diversity of the bacterial and archaeal community present in 
hummock and hollow soils of three geographically separate Sedge- 
dominated wetlands on the Qinghai-Tibetan Plateau. 454- 
pyrosequencing of 16S rRNA genes was used to analyze the 
relative abundance of major bacterial and archaeal groups. We 
also analyzed the correlation of environmental factors to the 
microbial community structures in this distinct ecosystem. 

Materials and Methods 

Sampling permissions 

Permission for sampling at the Hongyuan wedand was obtained 
from the Administrative Office of the Hongyuan Riganqiao 
Wetland Natural Reserves. No specific permissions were required 
for the sampling in Dangxiong and Maduo wetlands. The field 
studies did not involve endangered or protected species. 

Site description and soil sampling 

Three Sedge-dominant wetlands located on the Qinghai- 
Tibetan Plateau (Figure SI) and characterized by hummock- 
hollow microtopography were chosen for the study. The 
Dangxiong (DX) wetland (30°28'1 l.OON, 91°03'42.02E, 4290 m 
a.s.l.) is located in DX county of Tibet nearby the southern edge of 
the Nyainqentanglha Mountains and in the basin of two tributaries 
of the Lhasa River. The main vegetation is Kobresia and Carex in 
hummocks and hollows, respectively. The average annual air 
temperature is 1.3°C and the average annual rainfall is 476.8 mm 
[20]. The second wetland is the Riganqiao/Hongyuan (HY) 
peadand (33°06'52.73N, 102°38' 37.39E, 3459 m a.s.l.) located in 
HY County (northwestern Sichuan Province), and is located in the 
northeastern edge of the Qinghai-Tibetan Plateau. The climate of 
HY is a continental high-plateau monsoon climate, with mean 
annual air temperature of 1.1 °C and average annual precipitation 
of 650 mm [21]. The vegetation cover is primarily Carex 
meyeriana and Carex muliensis [22]. Both DX and HY wedands 
have mature peat habitats [21,23]. The third wetland is Maduo 
(MD) wedand (34°38' 10.30N, 98°01 '46.29E, 4229 m a.s.l.) 
located in the center of Sanjiang River Source Natural Reserve 
[7]. Typical peat has not been formed in MD wetland. The 
average annual air temperature of MD is — 4. 1°C and the mean 
annual precipitation is 303.9 mm [24]. The dominant plants in 
MD wedand are Kobresia (in hummocks) and Carex (in hollows). 



At the time of sampling, the water table was ca. 1 5 cm above 
hollow soil surfaces and ca. 5 cm below the top hummock soils. 
Soil cores from three hummocks and three hollows in these three 
wetlands were collected from 0-5 cm soil depth in August 201 1. 
Soil samples were kept in a cool box during transportation and 
were stored in the laboratory at 4°C. Soil pH was determined after 
transportation of samples to a laboratory using a compound 
electrode and a soil-to-water ratio of 1:2.5. Soil moisture (SM) was 
measured by drying soils at 105°C for 24 h. Soil organic carbon 
(SOC) and total N (TN) were determined by dichromate oxidation 
and Kjeldahl digestion, respectively. Available P (AP) in the soil 
was measured by the sodium bicarbonate extraction-molybde- 
num-antimony anti-colorimetry method. 

Some abbreviations are used in the following paragraphs and 
figures for simplicity: DXa and DXb represent the hummock soils 
and the hollow soils of the DX wedand, respectively; HYa and 
HYb represent the hummock soils and the hollow soils of the HY 
wetland and MDa and MDb represent the hummock soils and the 
hollow soils of the MD wedand. 

DNA extraction, real-time PCR and pyrosequencing 

Soil DNA was extracted from 0.4 g fresh soil (18 samples, 3 
hummock soils and 3 hollow soils from 3 wetlands) using the 
NucleoSpin soil kit (Macherey-Nagel, Diiren, Germany). DNA 
quality was verified using a NanoDrop 1000 spectrophotometer 
(Thermo Fisher Scientific, Schwerte, Germany) and diluted 100- 
fold in water for quantitative real-time PCR (qPCR) analysis. The 
assays targeting bacterial 16S rRNA and archaeal 16S rRNA were 
performed using the Taqman real-time PCR System as described 
previously [25]. The assays were performed using an iCycler 
instrument (Bio-Rad) and the associated software. DNA extracts 
were diluted 10-fold before the PCRs for pyrosequencing. Primers 
F515 (5 ' -GTGCC AGCMGC C GCGGTAA) and R806 (5'- 
GGACTACVSGGGTATCTAAT) were used to amplify 16S 
rRNA genes. Individual PCRs were barcoded by 6-bp molecular 
barcodes specific for each sample for the subsequent identification. 
PCRs were performed in 50-jLll volumes containing 5 |ll 
lOxAccuPrime PCR Buffer II (Life Technologies, Darmstadt, 
Germany), 0.4 mM of each primer, 1 (Lll of Taq AccuPrime (Life 
Technologies), 1 |il of template and 39 |ll sterile water. Cycling 
was performed with an initial denaturation at 94°C for 5 min 
followed by 27 cycles: 50°C 30 s, 68°C 30 s, 94°C 30 s, and a final 
extension at 68°C for 10 min on an Eppendorf Mastercycler 
instrument (Eppendorf, Hamburg, Germany). PCR products from 
each tagged primer set were pooled and purified using the 
GenElute PCR Clean-Up Kit (Sigma, Taufkirchen, Germany), 
and DNA concentrations were determined using a Qubit 
instrument (Life Technologies). Finally, samples were pooled in 
an equimolar concentration for 454-pyrosequencing. Pyrosequen- 
cing was carried-out at the Max Planck-Genome-Centre Cologne 
and performed using standard procedures using a Roche 454 
Genome Sequencer GS FLX+ instrument. The 454 pyrosequen- 
cing reads (raw data) were deposited under the study accession 
number SRP033622 in the NCBI Sequence Read Archive. 

Post-run sequence analysis 

The taxonomic assignment of pyrosequencing reads and 
preprocessing of sequences was performed using USEARCH (v. 
7.0.1090) [26], Qiime (v. 1.7.0) [27] and Mothur (v. 1.27) [28] 
software platforms, as described below. All raw sequences were 
first analyzed and the OTU table was created using the UPARSE 
pipeline [26] (http://www.drive5.com/usearch/manual/ 
uparse_cmds.html). In this pipeline, sequences were first sorted 
based on barcodes and sequences were removed from further 
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analysis when they were shorter than 200 bp, contained ambig- 
uous bases or homopolymers greater than 6 bp in length. 
Chimeras were removed using UCHIME [29]. Operational 
taxonomic units (OTUs) were picked from the good sequences 
with a cutoff value of 97% sequence identity. Taxonomic 
classification was carried-out in Qiime with the Bayesian classifier 
and the SILVA reference database with a confidence threshold of 
80%. Rarefaction curves and diversity indices including microbial 
community diversity (Shannon index), richness (Chaol) and 
coverage were calculated in Mothur. 

A heatmap representation of the relative abundance of OTUs 
between samples was constructed using R (http://www.r-project. 
org/). A Hellinger transformation of the OTU counts was 
performed using the decostand function in the vegan package 
[30]. Principal components analysis (PC A) was performed using 
prcomp and the result indicated that PCI, PC 2 and PC 3 
explained 25%, 15% and 12% of the variance, respectively. In 
order to select the OTUs explaining most of the differences 
between samples, the 25 OTUs with highest loadings of PCI, 
15 OTUs of PC 2 and 12 OTUs of PC 3 were chosen to construct 
the heatmap [31]. A total of 42 unique OTUs were obtained since 
some of the OTUs were selected from more than one PC. The 
OTU abundances were converted to percentage of reads from 
each sample and Manhattan distances were calculated and the 
heatmap constructed using the heatmap. 2 function in gplots [32]. 
The taxonomy of the selected OTUs was added separately and 
ANOVA was performed using SPSS Statistics (20.0) Software 
(IBM). 

UniFrac is a P-diversity measure that uses phylogenetic 
information to compare environmental samples. Unweighted 
UniFrac distances [33] in principle coordinate analysis (PCoA) 
were used to compare the microbial community composition 
among samples. Datasets were subsampled to 4,130 sequences. 
Representative OTU sequences from UPARSE were aligned in 
Mothur, and then a newick formatted phylogenetic tree was built 
using FastTree in Qiime. Lastly, a 3D PcoA plot was created in 
Qiime and visualized by KiNG (v. 2.16). 

Two-way analysis of variance (ANOVA) on microbial diversity 
indices and relative sequence abundances of major taxa and 
selected OTUs, and the subsequent statistical analysis of simple 
main effects and the post hoc pairwise test were all performed by 
SPSS statistics (20.0) software. To identify environmental factors 
influencing microbial diversity, canonical correspondence analysis 
(CCA) was performed between the OTU profiles (subsampled to 
4,130 sequences per sample) and environmental factors using the 
vegan package in R [30]. 

Results 

Microbial diversity indices 

A total of 158,939 microbial sequences belonging to 
9,907 OTUs were obtained after the quality filtering. A similarity 
level of 97% was used to identify OTUs and to estimate diversity. 
At this genetic distance, an average of 1,179 OTUs and Chaol 
value of 2,156 were found in the samples at a sequencing depth of 
4,130 (Table SI in File SI). The rarefaction curves of OTUs and 
Chaol were not saturated (Figure 1A and B). In contrast, 
rarefaction curves of Shannon indices approached a plateau at a 
sequencing depth of 2,000 (Figure 1C). The Shannon diversity 
indices at a sequencing depth of 4, 1 30 and its rarefaction curves 
showed that it ranged from 5 to 7 (Figure 1 C and Table S 1 in File 
SI). The coverage of all samples at a sequencing depth of 4, 130 is 
above 80% (Figure ID and Table SI in File SI). 



Two-way ANOVA was used to analyze the effects of locations 
and microtopographies on all diversity indices, and significant 
interaction effects were detected (Table SI in File SI, P<0.05). 
The simple main effects analysis showed that the main reason for 
these interaction effects was the low diversity and richness of DXb 
samples. OTUs, Chaol, Shannon and Coverage indices between 
hummock and hollow soils were only significandy different in the 
DX wetland (results not shown). 

The abundance of archaeal taxa 

The abundances of Archaea in these three wetiands were 
significandy different. HY wetland had the highest abundance of 
Archaea at about 10 16S rRNA gene copies g soil, whereas it 
was 10 6 -10 8 copies g" 1 soil in DX and only 10 6 -10 7 copies g" 1 
soil in MD wetland (Table 1). The hollow soils in DX and MD 
wetlands had significantly more archaeal sequences compared 
with the hummocks. 

The pyrosequencing PCR primers used for the amplification of 
16S rRNA genes have specificity for both bacterial and archaeal 
sequences. A total of 65 archaeal 16S rRNA OTUs from 2,549 
sequences were obtained, corresponding to 2.61%, 4.74%, 0.14%, 
0.96%, 0.07% and 0.49% of microbial sequences from HYa, 
HYb, DXa, DXb, MDa and MDb, respectively. These ratios were 
similar to the ratio calculated using the qPCR data (Table 1), 
which were 2.74%, 3.89%, 0.11%, 0.75%, 0.2% and 0.66%. 

Most samples of DXa and MDa had less than 1 0 sequences, so 
these were not included in the taxonomy analysis because the 
coverage was too low. The archaeal communities were composed 
of both Crenarchaeota and Euryarchaeota phyla (Figure 2). 
Crenarchaeota-affiliated sequences were relatively abundant in 
MDb (30.7%). The crenarchaeotal sequences were mainly 
affiliated with Miscellaneous Crenarchaeotic Group and Soil 
Crenarchaeotic Group. Euryarchaeota sequences were dominant 
in most of the samples and accounted for 97.8%, 93.7%, 84.7% 
and 69.3% of archaeal sequences from HYa, HYb, DXb and 
MDb, respectively. 

Relative abundance of various bacterial taxa 

The 16S rRNA gene sequences retrieved from all samples were 
affiliated with 36 bacterial phyla and corresponded to 9,772 
bacterial species-level OTUs. The relative abundance (>1%) of 
various phyla and proteobacterial classes is shown in Figure 3. 
The predominant phylum was Proteobacteria (average 34.5%), 
including Alphaproteobacteria (16.9%), Betaproteobacteria 
(7.2%), Deltaproteobacteria (5.8%) and Gammaproteobacteria 
(3.7%) classes. The community structures were significandy 
affected by location (P<0.001, Table S2 in File SI). MD soils 
had a significandy (P<0.001) higher relative abundance of 
Alphaproteobacteria sequences than DX and HY soils (Figure 3). 
Also, significantly higher relative abundances of Alphaproteobac- 
teria were obtained from hummocks than hollows (P = 0.004, 
Figure 3). Orders Rhizobiales and Rhodospirillales comprised 
78.6% of sequences in the class Alphaproteobacteria, and 
Rhodospirillales in particular had higher relative abundance in 
MD soils compared with the other two wetlands (Figure S2A). The 
most abundant phylotype within the Betaproteobacteria belonged 
to the Burkholderiales and the phylotype with the second highest 
relative abundance could not be classified at a higher level and was 
named "Others" in Figure S2B. Within Deltaproteobacteria, 
Myxococcales and Desulfuromonadales were the two orders with 
highest relative abundance (Figure S2C). Gammaproteobacteria 
sequences, especially Xanthomonadales, the most abundant order 
in this class, had significantly lower relative abundances in HY 
soils compared with the other two wetlands (Figure S2D). 
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Figure 1. Rarefaction curves showing the relationship between the increase of bacterial ((-diversity and the number of randomly 
sampled sequences. A, number of OTUs observed experimentally; B, diversity predicted by Shannon index; C, richness predicted by Chaol; D, the 
sample coverage. 
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The other phyla with high relative abundances, as shown in 
Figure 3, were Actinobacteria (17.3%), Bacteroidetes (11%), 
Chloroflexi (6.4%), Acidobacteria (5.5%), Verrucomicrobia 
(4.3%), Firmicutes (4.1%), and Planctomycetes (4%). Their relative 
abundances were significantly different between locations (P< 
0.05) except for Chloroflexi and Planctomycetes (Table S2 in File 
SI); however, the relative abundance of Planctomycetes was 
significantly higher in hummocks than hollows across all samples 
(P = 0.004, Figure 3) and Chloroflexi was significantly higher in 
hollows than hummocks in HY and MD samples (P = 0.005, 
Figure 3). In addition to the above abundant phyla, rare bacterial 
phyla were also found and are shown in Figure S3. We defined 
rare phyla as those having an average abundance between all 
samples of ^1%. Although Candidate_division_OP8 and Cyano- 
bacteria were rare when considering all samples together, their 
relative abundances were > 1 % in HYb and DXb respectively. In 
total, these rare phyla together constituted a higher proportion of 
the sequences in hollows compared with hummocks (Figure S3). 

Acidobacteria, Bacteroidetes, Chloroflexi and Verrucomicrobia 
phyla were relatively abundant (>27%) in the study sites. Within 
the Acidobacteria phylum, Acidobacteria and Holophagae were 



the two most abundant classes. The relative abundance of the 
Acidobacteria class was affected by location (P = 0.001, Table S3 
in File SI) and had low abundance in HY soils (Figure S4A). In the 
phylum Bacteroidetes, the most abundant classes were Sphingo- 
bacteria and Bacteroidia, which together amounted to 86.6% of 
the Bacteroidetes phylum (Figure S4B). Sphingobacteria had 
higher relative abundance in hollow soils compared with 
hummock soils (P= 0.004, Table S3 in File SI, Figure S4B). 
The two most abundant classes in the phylum Chloroflexi were 
Anaerolineae and KD4-96, contributing 75.1% to this phylum 
(Figure S4C). Anaerolineae was also associated with location 
(P = 0.001) and had greater relative abundance in hollow soils 
compared with hummock soils (P<0.001, Table S3 in File SI, 
Figure S4C). Within the Verrucomicrobia phylum, the classes 
Spartobacteria and OPB35 had the highest relative abundance 
(Figure S4D). Spartobacteria was significandy more abundant in 
hummocks than hollows across all sites (P<0.001, Table S3 in File 
SI, Figure S4D). 
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Figure 2. Taxonomic assignment of sequences affiliated with 
Archaea from samples in different wetlands and microtopo- 
graphies. Error bars indicate the standard deviation of relative 
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Methane cycle-related microorganisms 

Methanogens and methanotrophs were present in the hum- 
mocks and hollows of three wetlands. The HY wetland had the 
highest abundance of both methanogens and methanotrophs 
(Table 2). The main classes of methanogens were Methanobac- 
teria and Methanomicrobia, which together accounted for 94.6%, 
87.1%, 66.4% and 57.4% of all archaeal reads in HYa, HYb, 
DXb and MDb, respectively (Figure 2). Methanobacterium was the 
most abundant genus of methanogen in HY wetlands. Methano- 
saeta were only present in DX and HY wetlands, and tended to 
have higher relative abundance in hollow soils of the HY wetiand. 
Methanosarcina had higher relative abundances in hollows of the 



■mil 

DXa DXb HYa HYb MDa MDb 



Bacteria-others 

Planctomycetes 

Firmicutes 

Verrucomicrobia 

Acidobacteria 

Chloroflexi 

Bacteroidetes 

Acti no bacteria 
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Deltaproteobacteria 
i Betaproteobacteria 
i Alphaproteobacteria 



Figure 3. Relative sequence abundances of bacterial phyla and 
proteobacterial classes. Major taxa detected with average relative 
sequence abundances >1% are displayed. Column "Bacteria-others" 
indicates combined relative sequence abundances of all the rare phyla, 
candidate divisions and of the taxonomically unclassified sequences; 
rare phyla were defined as having average relative sequence 
abundances between all samples of <1%. 
doi:10.1371/journal.pone.0103115.g003 
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three wetlands, and also in hummocks of HY wetland. Methano- 
trophs detected in these wetlands belonged to the phylum 
Proteobacteria, including type I methanotrophs (Gammaproteo- 
bacteria, Methylobacter) and type II methanotrophs (Alphaproteo- 
bacteria, Methylocystis and Methylosinus). Methylocystis was 
dominant and had high relative abundance in HY soils compared 
with the other two sites. 

OTU-level microbial p-diversity analysis 

P-diversity assesses the differences between microbial commu- 
nities and reflects the dissimilarity between samples. A heatmap 
was used to intuitively display the differences in OTU relative 
abundances observed in this study (Figure 4). The OTUs with 
highest contribution to the PCA ordination were selected. Most of 
these OTUs (36/ 42) were significantly different between the three 
locations based on ANOVA tests (Figure 4). The OTUs fell into 
five main clusters. The cluster positioned at the top of the heatmap 
mostly had high relative abundance in the MD wetland only. 
Among them, eight OTUs belonged to the class Alphaproteo- 
bacteria (five OTUs of Rhodospirillaceae and three of Sphingo- 
monadaceae). The three clusters of OTUs in the middle of the 
heatmap were mostiy associated with high relative abundance in 
the HY and DX wetlands. Within this group was OTU- 17, 
corresponding to the genus Methanobacterium. The cluster at the 
bottom of the heatmap comprised OTUs that had high relative 
abundance in the hollows of the DX wedand. The relative 
abundance of most OTUs (36 of 42) exhibited significant 
differences based on location, and many (26 of 42) based on 
topogeography. For example, OTU- 7 (belonging to Bradyrhizo- 
biaceae), OTU-1 (Micrococcaceae) and OTU-4 (Sphingobacter- 
iales, vadinHA17) had higher relative abundance in hollows 
compared with hummock soils (P<0.01). 

A PCoA plot was used to visualize the data based on P-diversity 
metrics of unweighted UniFrac (Figure 5). There were obvious 
separate, non-overlapping groups for each of the locations, and 
hummock samples were also separated from hollow samples. A 
partial Mantel test was used to estimate the correlation between 
the OTU level unweighted UniFrac matrix and another matrix 
derived from either the location or microtopography parameters, 



while using one of the other location or microtopography as the 
third "control" distance matrix. The Mantel r statistic of 0.55 and 
P-value of 0.001 indicated a significant relationship between 
UniFrac distances and wetland locations while controlling for 
differences in microtopographies. There was also a relatively 
strong relationship between UniFrac distances and microtopo- 
graphy while controlling for differences in locations (r = 0.41, 
P = 0.001). These results confirmed that bacterial communities 
were significantly different between locations and microtopogra- 
phies. 

Association between microbial community structures 
and soil variables 

The environmental data of the three wetlands are shown in 
Table 1. We used Mantel's tests to determine the significance of 
independent influences of environmental factors (AP, SOC, TN, 
C:N, pH, SM, plant biomass (PB), mean annual air temperature 
(T), mean annual precipitation (P) and altitude (A)) on the 
distribution and diversity of microbial OTUs across different 
samples. Except plant biomass, significant correlations were found 
between environmental factors and microbial community struc- 
ture. The soil factors (AP, SOC, TN, pH and SM) and the location 
factors (T, P and A) all explained > 1 0 % of microbial community 
variation between samples (Table S4 in File SI). Canonical 
correspondence analysis (CCA) was used to discern the microbial 
structure according to environmental parameters. The seven 
significant environmental variables were selected and included in 
the CCA plot (Figure 6). pH was the most important contributor 
to axis 1, which on its own accounted for 13.7% of the variation. 
These environmental variables together explained 64% of the 
difference in the microbial communities between the sites. 

Discussion 

We analyzed the microbial diversity in different wetlands of the 
Qinghai-Tibetan Plateau using identical experimental conditions, 
which is critical for taxon richness estimation. The microbial 
rarefaction curves of OTUs and Chaol were not saturated 
(Figure 1A and B), suggesting we still did not survey the full extent 



Table 2. Percent relative abundance of 16S rRNA genes in the pyrosequencing datasets associated with genera of methanogens 
and methanotrophs (mean±SD, n = 3). 





Taxon 


DXa 


DXb 


HYa 


HYb 


MDa 


MDb 


Methanogens 


Methanobacterium 


0.05±0.07 


0.29±0.18 


1.62 ±0.21 


2.56±1.51 


0.02 ±0.01 


0.04±0.03 


Methanocella 


0 


0.02±0.01 


0.02 ±0.02 


0.05 ±0.04 


0 


0 


Candidatus Methanoregula 


0±0.01 


0±0.01 


0.07±0.08 


0.27±0.18 


0 


0 


Rice Cluster II 


0 


0 


0.05 ±0.04 


0.11 ±0.05 


0 


0 


Methanosaeta 


0.01 ±0.01 


0.06±0.09 


0.06±0.02 


0.76±0.37 


0 


0 


Methanosarcina 


0.01 ±0.01 


0.26±0.23 


0.60±0.54 


0.25±0.12 


0 


0.21 ±0.16 


Methanotrophs 


Methylocystis 


0.19±0.32 


0.09 ±0.06 


0.69±0.15 


0.62±0.12 


0.03 ±0.03 


0.04±0.05 


Methylosinus 


0 


0±0.01 


0.02±0.01 


0±0.01 


0 


0.01 ±0.02 


Methylobacter 


0 


0±0.01 


0.02±0.02 


0.07±0.08 


0 


0.01 ±0.01 


Methylocaldum 


0 


0 


0.03±0.04 


0.05±0.01 


0 


0 


Methylomonas 


0.01 ±0.01 


0 


0.02±0.01 


0.02±0.02 


0 


0 



Only genera with >0.01% relative abundance are shown. 
doi:1 0.1 371 /joumal.pone.01 031 1 5.t002 
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Figure 4. Heatmap showing the relative abundance of selected OTUs. The samples and OTUs are clustered according to Manhattan 
distances. The colors correspond to the relative abundance of the OTUs in the samples, as indicated by the color legend. The taxonomy of each OTU 
is provided to the lowest-level attained during the classification. Abbreviations are used to indicate class (c), order (o), suborder (so), family (f) and 
genus (g). Two-way ANOVA for each OTU was performed using location (L), topography (T) and location and topography (LxT). P-values of the 
ANOVA are indicated as follows: * (P<0.05), ** (P<0.01), *** (P<0.001). 
doi:1 0.1 371 /journal.pone.01 031 1 5.g004 



of taxonomic diversity at a taxonomic resolution of 97% sequence 
identity for the fragments of 16S rRNA genes analyzed. At the 
sequencing depth of 4,130, the Chaol and average detected 
OTUs were lower in hollow soils of DX wetland than other 
samples; however, even these indices in DXb (Chaol value of 
1,563 and 8,75 OTUs, Table SI in File SI) were higher than those 
reported for a northern peatland (Chaol value of ~700 and 
~500 OTUs) [2]. The rarefaction curves of the Shannon index 
approached the plateau phase with more than 2,000 sequences per 
sample, indicating that Shannon curves were saturated (Fig- 
ure IB). The differences in Shannon indices between the three 
wetlands were not significant; however, relatively lower Shannon 
indices (6.2) were found in our study site compared with other 
natural wedands, e.g. 7.3 in northern peatland of Russia [2] and 
6.7 in Yellow River Estuary wetlands [34]. Although we analyzed 
the sequences with quality controls, some errors may still exist 
because of PCR error, sequencing error, and the methods 
employed within the bioinformatic pipeline [35], Although the 
Shannon index is relatively insensitive to sequencing and PCR 
errors [35], it is still sensitive to the bioinformatic processing. 
Therefore, even with the same high-throughput sequencing 
method, it is difficult to directly compare the diversity indices 
obtained using different protocols. 

The large bacterial diversity in these wetlands included high 
relative abundances of Proteobacteria (37.5%), Actinobacteria 
(17.3%) and Bacteroidetes (11%) (Figure 3). It is not surprising to 
see the dominance of Proteobacteria and Actinobacteria as they 
are prevalent in soil of various ecosystems [36]; however, 
Bacteroidetes had the third highest relative abundance among 
the bacterial phyla and these were greater than those of wetlands 
in most other regions. For example, wedands of Northern Russia 
[2], a North Carolina coastal plain and the Florida Everglades 
[19] all had Bacteroidetes totaling less than ~1%. Bacteroidetes 
were found to have the highest relative abundance among 
microorganisms in high-altitude wetlands in Chile [37]. Also, 
they were found to be the predominant bacteria in Tibetan 

PC 2 (11%) 



^IDa2jtfDai 



.MDbl DX&ljpXb2 

.DXb 3 
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HYa ,/IYa2 JHYb2 
^HYb3 
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plateau glaciers [38], and also one of the dominant phyla in 
sediments [17] and freshwaters [39,40,41] of the Tibetan plateau 
lakes. Therefore, the presence of relatively high abundance of 
Bacteroidetes in the three wetlands of the Qmghai-Tibetan 
Plateau may be partly associated with factors shared between 
these high altitude environments and warrants further investiga- 
tion into their ecological role in these systems. 

The relative abundance of Acidobacteria was low in the 
Qinghai-Tibetan wetland samples compared with some other 
wedands. For example, Acidobacteria was the dominant phylum 
in a northern peatland [2] and a bog of the glacial Lake Agassiz 
peatland [4]. These two sites have pH values of ~4. In contrast, 
Acidobacteria was not highly abundant in the fen sediments of the 
Lake Agassiz peadand [4], which had a slightly higher pore water 
pH value of ~6. The average porewater pH of the HY and DX 
peadand was 6.3 and in MD it was 7.8, which could explain the 
relatively low Acidobacteria abundances observed. In agreement 
with our results, a previous study also did not detect Acidobacteria 
by DGGE in peatland of Zoige plateau close to where HY wedand 
is located [14]. 

Methane in natural wedands is the final product of the 
anaerobic degradation of organic matter and performed by 
methanogenic archaea [42] . A large part of methane is oxidized 
by methanotrophs before it reaches the atmosphere [43] and 
therefore both methanogens and methanotrophs are the key 
players in the methane-cycling of natural wetlands. Methanogens 
in HY wetland were significantly more abundant than the other 
two sites, which is consistent with the higher SOC content in this 
wedand. Methanogens were more abundant in the hollows, which 
is consistent with their strict anaerobic metabolism and the hollows 
more likely to be anoxic compared with hummocks. Methano- 
bacteria were the most abundant methanogens in HY and DX 
wetlands, followed by Methanomicrobia. The opposite was found 
in the MD wetland with higher relative abundance of Methano- 
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Figure 5. Principal coordinate analysis (PcoA) plot of samples 
using the unweighted UniFrac distance metric. 

doi:1 0.1 371 /journal.pone.01 031 1 5.g005 



CCA1 (13.9%) 

Figure 6. Canonical correspondence analysis (CCA) used to 
evaluate the effect of soil factors on the bacterial community 
structure. Soil factors indicated in blue text include AP (available 
phosphorus), pH, SOC (soil organic carbon), TN (total nitrogen), T (mean 
annual air temperature), SM (soil moisture) and A (altitude). 
doi:10.1 371/journal.pone.01 031 1 5.g006 
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microbia followed by Methanobacteria. Previous reports from the 
Zoige wedand on the Qinghai-Tibetan plateau identified members 
of Methanomicrobia (Methanosarcinales and Methanomicrobiales) 
as the dominant methanogens [8,44]. This variability suggests that 
the methanogen community composition in these wetlands might 
be a result of multiple factors and further studies should be 
performed to determine which factors play key roles in structuring 
methanogen populations. Among methanotrophs, Methylocystis 
was the most abundant, which is consistent with the results of a 
pmoA gene pyrosequencing study of the HY site [10]. The pmoA 
gene encodes a subunit of the methane monooxygenase enzyme 
and is frequently used to identify methanotrophs in environmental 
samples. A comparison between these studies showed relatively 
fewer Type I methanotrophs (Gammaproteobacteria) detected by 
16S rRNA than pmoA gene pyrosequencing. This can be 
explained by many of the Type I-associated pmoA clades being 
associated with uncultivated methanotrophs [10], which in many 
cases would not be identified based on 16S rRNA sequences 
depending on their phylogenetic distance to known methano- 
trophs. 

Our analysis of the relative abundance of various microbial 
phyla, classes, orders and especially OTUs, showed that the major 
microbial community structures were significandy different among 
wetland locations and microtopographies. The data indicated that 
pH was one of the major soil properties influencing the microbial 
community composition and diversity in wetlands of the Qinghai- 
Tibetan Plateau. This effect of environmental pH has also been 
found in other terrestrial and aquatic ecosystems [16,17,45]. The 
influence of pH on microbial community composition is probably 
due to the narrow pH ranges for optimal growth of bacteria [16]. 
Among the soil properties, SOC, AP, TN, and SM were also 
important contributors to the variations in the microbial 
community within these wedand soils. All of these variables are 
related to nutrient availability, which will have obvious implica- 
tions for microbial and plant growth. A total of 64% of community 
variance could be explained by these environmental factors, 
leaving 36% unexplained. Therefore, additional factors such as 
plant species [46] and some unmonitored variables are likely to 
have further contributed to the differences in microbial commu- 
nity structure. 

Conclusions 

In summary, microbial communities of three Sedge-dominated 
natural wetlands in the Qinghai-Tibetan Plateau were investigated 
by 454-pyrosequencing in this study. The differences of soil 
microbial community structure between the wetiands at different 
geographic locations and microtopographies could be clearly 
distinguished by pyrosequencing. The CCA results suggested that 
there was a clear association between microbial community 
structure and wedand soil properties. These results provide insight 
into the microbial community structure in this distinct ecosystem 
and identified the main factors shaping microbial community 
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